An Axisymmetric Gravitational Collapse Code 



Matthew W. Choptuik 
CIAR Cosmology and Gravity Program, Department of Physics and Astronomy, 
University of British Columbia, Vancouver BC, V6T 1Z1 Canada 

Eric W. Hirschmann 

Department of Physics and Astronomy, Brigham Young University, Provo, UT 84604 

Steven L. Liebling 
Southampton College, Long Island University, Southampton, NY 11968 

Frans Pretorius 

Theoretical Astrophysics, California Institute of Technology, Pasadena, CA 91125 

£C) , We present a new numerical code designed to solve the Einstein field equations for axisymmetric 

f*"^ ■ spacetimes. The long term goal of this project is to construct a code that will be capable of studying 

f*"*) ' many problems of interest in axisymmetry, including gravitational collapse, critical phenomena, 

investigations of cosmic censorship, and head-on black hole collisions. Our objective here is to detail 
the (2+l) + l formalism we use to arrive at the corresponding system of equations and the numerical 
methods we use to solve them. We are able to obtain stable evolution, despite the singular nature 
of the coordinate system on the axis, by enforcing appropriate regularity conditions on all variables 
and by adding numerical dissipation to hyperbolic equations. 
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J> . I. INTRODUCTION 

. 

In this paper we introduce a numerical code designed to solve the Einstein field equations for axisymmetric 
spacetimes. Even though the predominant focus in numerical relativity in recent years has been to study 
situations of relevance to gravitational wave detection, and hence lacking symmetries, there are still numerous 
ff) ' interesting problems, both physical and computational, that can be tackled with an axisymmetric code. The 

advantages of restricting the full 3D problem to axisymmetry (2D) are that the complexity and number of 
equations are reduced, as are the computational requirements compared to solving a similar problem in 3D. 

Prior numerical studies of axisymmetric spacetimes include head-on black hole collisions 0, S EL El i 
collapse of rotating fluid stars 0, IS U EJ > the evolution of collisionless particles applied to study the stability 
of star clusters [llj and the validity of cosmic censorsh ip Il 2 | . evolution of gravitational waves 
black hole-matter-gravitational wave interactions |15L iWL \YA Il8| , and the formation of black holes through 
gravitational wave collapse an d corresponding critical behavior at the threshold of formation [2(j • 

Our goals for creating a new axisymmetric code are not only to explore a wider range of phenomena than 
those studied before, but also to provide a framework to add adaptive mesh refinement (AMR) and black 
hole excision to allow more thorough and detailed investigations than prior works. 

The outline of the rest of the paper is as follows. In Section|n]we describe the (2 + 1) + 1 decomposition 0, 
22^ of spacetime that we adopt to arrive at our system of equations. The (2 + 1) + 1 formalism is the familiar 
ADM space + time decomposition (in this case 2+1) applied to a dimensionally reduced spacetime obtained by 
dividing out the axial Killing vector, following a method devised by Geroch 23]. In Section ITTT1 we specialize 
the equations to our chosen coordinate system, namely cylindrical coordinates with a conformally flat 2— 
metric. At this stage we do not model spacetimes with angular momentum, and we include a massless scalar 
field for the matter source. In Section llVl we discuss how we search for apparent horizons during evolution. In 
Section we describe our numerical implementation of the set of equations derived in Section IlIII A variety 
of tests of our code are presented in Section fVTl which is followed by conclusions in Section IvTTl Some details 
concerning our finite difference approximations, solution of elliptic equations via the multi-grid technique, and 
a spherically symmetric code used for testing purposes are given in Appendices A and B. Unless otherwise 
specfied, we use the units and conventions adopted by Misner, Thorne and Wheeler [24|. 
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II. THE (2+l)+l FORMALISM 



The most common approach in numerical relativity is to perform the so-called 3 + 1, or ADM, split of the 
spacetime that one would like to evolve. In this procedure, a timelike vector field is chosen together with 
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spatial hypersurfaces that foliate the spacetime. If axisymmetry is assumed, it is usually incorporated once 
the ADM decomposition has been done, and is reflected in the independence of the various metric and matter 
quantities on the ignorable, angular coordinate, (p. 

Our approach, which follows Geroch [2^| and Nakamura et al 0, reverses this procedure. We assume 
axisymmetry from the outset and perform a reduction of the spacetime based on this assumption. Once 
we have projected out the symmetry, we perform an ADM-like split (now a 2 + 1 split) of the remaining 
3-manifold. 

More specifically, we begin with a 4-dimensional spacetime metric 1 7 M „ on our manifold A4. The axisym- 
metry is realized in the existence of a spacelike Killing vector 

**-(*)'• 

with closed orbits. We define the projection operator g^, allowing us to project tensors from the 4- 
dimcnsional manifold yVf with metric 7^ to a 3-dimensional manifold A4/S 1 with metric g ab , as 



where s is the norm of the Killing vector 



With the definition of the vector 



X»X„ 



the metric on the full, 4-dimensional spacetime can be written as 



9ab + s z Y a Y b s'Y a 



s 2 Y b 



(2) 
(3) 
(4) 

(5) 



Projecting, or dividing out the symmetry amounts to expressing 4-dimensional quantities in terms of quan- 
tities on the 3-manifold. For instance, the connection coefficients ^r^„ associated with the 4-dimensional 
metric are 

(4) r A _ (3) r A , Q A 

= (3) ^, + \s 2 g X ° [YuZ^ + Y„Z va - d a (In s 1 ) Y„Y V ] + \y a [0„ (s 2 Y u ) + d v , (6) 

where ^r^„ are the connection coefficients constructed from the 3-metric g M „, and we have defined the 
antisymmetric tensor 



Zpu = %Y V - d v Y^. (7) 

Notice that is an intrinsically 3-dimensional object, in that Z^X^ = Z^X" = 0. With some algebra, 
the Ricci tensor on the 4-manifold, ^Raui can now be written as the Ricci tensor on the reduced space, 
^Rfiu, together with additional terms involving fields coming from the dimensional reduction 

wr^ = + d x q* v - + n^n^ - (8) 

where is the covariant derivative on the 3 dimensional manifold. Expressed in terms of s, Y a , Z ab and 
^R a b, the components of JSJ) are 
1 



1 L <pa 



-s^Z hc Z bc - sD a D a s, 



1 



2s 



D c {s A Z ac )+Y a 



s 4 Z bc Z bc - sD a D a s 



D n D h s 



1 



-S ZarZh 



■D c [ S 3 Z c{a \Y K 



Y,Y h 



-s"Z hc Z hc - sD a D a s 



(9) 
(10) 

(11) 



1 Note that we will use Greek letters to denote coordinates in the full 4— manifold and Latin letters to denote coordinates in the 
reduced 3-manifold. 
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Taking the trace of J8J , and using the definitions described above, gives the decomposition of the Ricci scalar 
as 

< 4) i? = <®R - -D a D a s - ]s 2 Z bc Z bc . (12) 
s 4 

The 4-dimensional Einstein equations, with a stress-energy tensor T^„, are 

(4) iV - \ {4) Rl^ = 87rT^. (13) 
Using equations 1)91121) . we can write the Einstein equations as 

8tt / 1 



D"D a s = -^-^{T w --sVy (14) 

D[ a w fc] = 8nse ab c T clp , (15) 

(3) i? afc = - s D a D b s + X -s 2 Z ac Z h c + 8ir \T^ a g v h - ^g ab T x x ^j , (16) 

where we have introduced the twist u> M of the Killing vector 

= S - e^Y" Z Xa , (17) 



and the four and three dimensional Levi-Civita symbols e^Acr and e a tc, respectively. The twist vector u> M is 
intrinsically 3-dimensional, i.e. w^X^ = 0. Furthermore, u> a /s 3 is divergence free 



tt'u 



= 0. (18) 



At this point, the first reduction is essentially done. Equation i|ltj[) can be viewed as the 3-dimensional 
Einstein equations, coupled to the projection of the stress-energy tensor T M „ and to induced "matter fields" 
s and w a (or Z ab ). The equations of motion for s and w a are given by 114JI and (jl5[l respectively; additional 
equations of motion will need to be specified for whatever true matter fields one incorporates into the system. 
This procedure so far is completely analogous to the Kaluza-Klein reduction from five to four dimensions, 
in which the 5-dimensional geometry becomes gravity in 4 dimensions coupled to electromagnetism and a 
scalar field |26j |. Here, however, the reduced 3-manifold has no dynamics and we have the rather appealing 
picture of having "divided out" the dynamics of the four dimensional gravitational system and reinterpreted 
the two degrees of freedom in the gravitational field as scalar (s) and "electromagnetic" (Z ab , or w a ) degrees 
of freedom. 2 

We now perform the ADM split of the remaining spacctime. This is done by first foliating the 3-dimensional 
spacetime into a series of spacelike hypersurfaces with unit, timelike normal vector n a . Then, similar to the 
dimensional reduction above, we decompose quantities into components orthogonal and tangent to n a , using 
the projection tensor h ab 

Kb = fjab + n a n b . (19) 

We now define the components of n a , and the induced 2-dimensional metric Hab 3 , using the following 
decomposition of the 3-metric 

g ab dx a dx b = -a 2 dt 2 + h AB (dx A + (3 A dt){dx B + (3 B dt), (20) 

where a is the lapse function and j3 A the shift vector. The gravitational equations now become the evolution 
equations for the components of the 2-metric Hab and the 2-extrinsic curvature Kab 

dthAB = -2aK A B + AApB + A B (3A (21) 
d t K AB = (3 c A c K A b + K A c^bP C + K BC A A P C +a[KK AB + (2) Rab] 

-2aK AC K B c - A A A B a-a^R A B, (22) 



2 Recall that the electromagnetic field in 3 dimensions has only a single degree of freedom as compared to the two degrees of 
freedom in 4 dimensions. 

3 We use upper-case Latin indicies to denote 2-dimensional tensor components. 
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the Hamiltonian constraint equation 

- K A B K B A + K 2 = ^R + 2 ^R ab n a n b (23) 

and the p and z momentum constraint equations 

A A K B A - A A K = - ^R cd n d h c B . (24) 

In the above, is the covariant derivative compatible with the 2-metric h A s, K = K A A , and ^'Rab 
and ( 2 >R are the 2-dimensional Ricci tensor and Ricci scalar, respectively. Note that because gravity in 3 
dimensions has no propagating degrees of freedom, the constraint equations fix the 3-dimensional geometry 
completely. Thus, if desired, one can use the constraint equations (|23I24|) instead of the evolution equations 
(I21I22|) to solve for h A B- The freely specifiable degrees of freedom of the 4-manifold are encoded in s and w a , 
which are evolved using 114|) . H15|) and (|18|l . Note that l|15fl and (|18() constitute four equations for the three 
components of w a — the purely spatial have four equations for the three components of w a - The purely spatial 
part of H15(l is, in 3 + 1 language, the angular momentum constraint equation and only needs to be solved at 
the initial time in a free evolution of w a . We note that the restricted class of axisymmetric spacetimes having 
no angular momentum (rotation) is characterized by the existence of a scalar twist, iu, such that w a = w^ a . 
In the vacuum case, w generally represents odd parity gravitational waves, while s encodes even parity, or 
Brill waves[? ]. We further note that this class includes the special case w a = 0, which will be the focus of 
our discussion below. 



III. COORDINATE SYSTEM, VARIABLES AND EQUATIONS 

In this section we describe a particular coordinate system and set of variables which, in the context of 
the formalism described in the previous section, provides us with the concrete system of partial differential 
equations that we solve numerically. We also detail the outer boundary conditions we use, and the on-axis 
regularity conditions necessary to obtain smooth solutions to these equations. 

We only consider spacetimes with zero angular momentum, and no odd-parity gravitational waves; therefore 
w a = 0. We choose a conformally flat, cylindrical coordinate system for the 2-metric 

h AB dx A dx B = tp(p, z, tf{dp 2 + dz 2 ). (25) 

This choice for h A s exhausts the coordinate freedom we have to arbitrarily specify the two components of 
the shift vector — (3 p (p,z,t) and f3 z (p,z,t). In order to maintain the form 125|) during an evolution, we use 
the momentum constraints, which are elliptic equations, to solve for (5 P and /3 Z at each time step. The 
Hamiltonian constraint provides a third elliptic equation that we can use to solve for the conformal factor i\>. 
For a slicing condition, we use maximal slicing of t = const, hypersurfaces in the ^-dimensional manifold — i.e. 
we impose = 0, where is the trace of the extrinsic curvature tensor of t — const, slices of 7 M „. This 
condition (specifically d^ 3 'K/dt = 0) gives us an elliptic equation for the lapse. 

Instead of directly evolving the norm of the Killing vector, s, we evolve the quantity a, defined by 

s = pi\?e pS 1 (26) 

and furthermore, we convert the resultant evolution equation for a (|14[) to one that is first order in time by 
defining the quantity SI, which is "conjugate" variable to a, via 

pCl = -2K/ - K z z 

= -^ Q ( lnS )- + ^^- ( 27 ) 

Part of the motivation behind using a and f2 as fundamental variables is to simplify the enforcement of on-axis 
regularity. In particular, regularity as p — * implies that a and O must exhibit leading order behavior of the 
form a = <J\{z, t)p + 0(p 3 ) and f2 = ^1(2;, t)p + 0(p 3 ) respectively, and experience has shown it to be easier 
to enforce such conditions, than to enforce the leading order behavior of s (or its time derivative) near the 
axis, which is s — s 2 (z, t)p 2 + 0(p 4 ). 

As mentioned previously, the only matter source we currently incorporate is a massless scalar field 3>(p, z, t), 
which satisfies the usual 4-dimensional wave equation □$ = 0. We convert this equation to first-order-in-time 
form by defining a conjugate variable II: 



IT = if) 2 n a $ 



(28) 
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The stress-energy tensor for the scalar field is 



(29) 



Using all of the above definitions and restrictions within the formalism detailed in the previous section, we 
end up with the following system of equations that we solve with our numerical code, described in the next 
section. The maximal slicing condition results in the following elliptic equation for a 



2 (pa, p ) p2 + a, zz + a p ^2-^ + {pa) 4 



4> A - 

2a 



•0 4 



6 a 



(8 P , P -B Z , Z ) 2 + (8 P , Z + B Z , P ) 2 
The Hamiltonian constraint gives an elliptic equation for tp 



2apCl + 8» iP - B z 



167raIT 



(30) 



%+8% + 16^ + 8(p«x) % + 8(p<r)^ 



1> 



1> 



1/, 



2a 2 



{B p , p -B z , z f + {B p , z + 8 z , p f 



6a 2 



2apCl + 8 P )P - B z 



-Wn (it 2 + $ p 2 + $ 2 2 ) 



6 (^(^),p),p3 

-2((pa), p ) 2 -2( p( T)^-2((^) J 2 . (31) 
The p and z momentum constraints, which provide elliptic equations that we use to solve for j3 p and /3 Z , are 



-0 P +0P + -8 Z 



ail 



■0 



(/3 p , P -/3%) 



6^f " (p*) 2 (/? P , 
a ip ' i 



(32) 



and 



-8 P 

3 P ' zp 



2ap 



%-*$-1<p*)*](p*-p*) + 



60— + 2 + 3fi(pa) 
2a / p0 6 



6 V a / ,p ; 



+ (pa) (f3 p , z + 8\ p ) = 



(33) 



The definition of il in Eq. (|27|l gives an evolution equation for a (where the overdot denotes partial differen- 
tiation with respect to t) 



a = 28" (pa)^+B z a r 



nil- 



P ' ,p 



(34) 



The evolution equation for il is 



n = 



2B p ( p n) p2+/ 3*n z - — (b 



z 2_ f3 p2 



4 \ p )-p 



0' 



v> 4 



V o ),p i/M 



4 V "0 



a.. 



(a. z 2tp z \ — 2 i 



647r^p($ !p2 ) 2 . 



(35) 



We also have an evolution equation for tp, which we optionally use instead of the Hamiltonian constraint l|31|) 
to update ip 



The definition of II and the wave equation for $ give 

$ = /5 p $ p + B z $ z 



2 



(36) 



(37) 



G 



and 



ft = pPU p + f3 z tt z + \u (apti + 2f3 p p + (3 Z Z ) 



1 

"ft 



2 (paft$ p ) p2 + (a^ 2 <P z ) z \ + ^ [(pa\ p $ p + (pa) z . (38) 



To complete the specification of our system of equations, we need to supply boundary conditions. In our 
cylindrical coordinate system, where p ranges from p = to p — p max and z ranges from z m - m to z max , we 
have two distinct boundaries: the physical outer boundary at p = p max , z — Zminj and z — z max ; and the axis, 
at p — 0. Historically, the axis presented a stability problem in axisymmetric codes. We solve this problem 
by enforcing regularity on the axis, and, as described in Section adding numerical dissipation to evolved 
fields. 

The regularity conditions can be obtained by inspection of the equations in the limit p — > 0, or more 
formally, by transforming to Cartesian coordinates and demanding that components of the metric and matter 
fields be regular and single valued throughout |25j . Garfinkle and Duncan 0] have further proposed that in 
order to ensure smoothness on the axis, one should use quantities that have either even or odd power series 
expansions in p as p — > 0, but which do not vanish faster than O(p). It is interesting that the quantities which 
we found to work best also obey this requirement. As discussed earlier, the particular choice of a and fl as 
fundamental variables was partly motivated by regularity concerns. The results are 

a, p (0,z,t) = (39) 

t P (0,z,t) = (40) 

P z , p (0,z,t) = (41) 

(3P(0,z,t) = (42) 

a(0,z,t) = (43) 

0(0, z,t) = (44) 

$ p (0,z,i) = (45) 

n iP (0,*,i) = (46) 

At the outer boundary, we enforce asymptotic flatness by requiring 

lim a(r,i) = l + ^ + 0(r- 2 ) (47) 
lim VM) = l + ^ + 0(r- 2 ) (48) 

r—*oc t 

lim p z (r,t) = ^ + 0(r- 2 ) (49) 

r—*oo f 

lim (3 p (r,t) = ffi + 0(r- 2 ), (50) 

r— >oo T 

for undetermined functions C(t), D(t), E(t), F(t), and r 2 = p 2 + z 2 . These latter relations are converted to 
mixed (Robin) boundary conditions (see Appendix A for details) and then are imposed at the outer boundaries 
of the computational domain: p = p max , z = z max , and z — — z max . We have also experimented with the use 
of Dirichlet conditions on a,(3 p and f3 z at the outer boundaries (specifically a = 1 and (3 P — (3 Z — there), 
and have found that these work about as well as the Robin conditions. For the scalar field, we assume that 
near the outer boundary we can approximate the field as purely radially outgoing, and require 

(r$) it + (r$) )f . = 0. (51) 

For scalar field configurations far from spherical symmetry, this approximation suffers and reflections are 
relatively large. However, in general, the reflections do not grow and are somewhat damped. For the other 
two evolved quantities, a and £1, we use this same naive condition for lack of any better, more physically 
motivated conditions. While this condition proves to be stable with damped reflections, a better condition is 
sought and this issue remains under investigation. 

For initial conditions, we are free to set (ct(0, p, z), 0(0, p, z), $(0, p, z), 11(0, p, z)). Once the free data is cho- 
sen, we then use the constraint and slicing equations to determine (a(0, p, z),ip(0, p, z), /3 Z (0, p, z), /3 P (0, p, z)). 
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Specifically, we define a general pulse shape 

\J {p- Pxf + e x (z - z x f - Rx 



Gx(p, z) = A x exp 



(52) 



characterized by six parameters (A x , px, ex, zx, Rx, Ax) and then choose initial data of the form 

5(0, p,z) = pG s {p,z) 

0(0, p,z) = pG n (p,z) 

$(0,p,z) = G$(p,z) 

H(0,p,z) = 0. (53) 

For ex = 1, these pulses are Gaussian, spherical shells centered at (px, zx) with radius Rx and pulse width 
Ax- For ex = 1 and (px, zx) = (0, 0), the pulses are spherical. The factor of p in the initial data for a and 
f2 ensures the correct behavior on axis for regularity. For the evolutions presented here, we let II = so that 
the initial configuration represents a moment of time symmetry. We note, however, that we are also able to 
generate and evolve approximately ingoing initial data. 

IV. FINDING APPARENT HORIZONS 

In this section we describe the equation and technique we use to search for apparent horizons (AHs) within 
t = const, spatial slices of the spacetime (see |2JJ |28| for descriptions of some of the methods available to 
find AHs in axisymmetry) . We restrict our search to isolated, simply connected AHs. In axisymmetry, such 
an AH can be described by a curve in the (p, z) plane, starting and ending on the axis at p = 0. We define 
the location of the AH as the level surface F = Q, where 

F = f-R(9), (54) 

and 



f = vV + (z- z ) 2 , (55) 
fsin^ = p, (56) 
f cos# = (z - zq). (57) 

The AH is the outermost, marginally trapped surface; hence, we want to find an equation for R(9) such that 
the outward null expansion normal to the surface F = 0, is zero. To this end, we first construct the unit 
spatial vector s a , normal to F = const. 



g ab F, 



(58) 



Then, using s a and the t = const, hypersurface normal vector n a , we construct future-pointing outgoing (+) 
and ingoing (— ) null vectors as 

l a ± = n a ± s a . (59) 

The normalization of the null vectors is (arbitrarily) t\t- a = —2. The outward null expansion 9+ is then the 
divergence of l\ projected onto F = const. 

9+ = (g ab - s a s b ) V b £ +a . (60) 

Using the definition of the extrinsic curvature Kab, and substituting l|59|l into l|60l) . we arrive at the familiar 
form for the null expansion when written in terms of ADM variables 

9+ = s A s A K AB + A A s A - K. (61) 

Note that because the normalization of 1°^ is arbitrary, so (to some extent) is that of 9 + . The above normal- 
ization is chosen so that 9 + measures the fractional rate of change of area with time measured by an observer 
moving along n a . 
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Substituting (|54|l and (|58|l into i|6f|) . and setting 9+ = 0, we are left with an ordinary differential equation 
for R(9). This equation takes the following form, where a prime ' denotes differentiation with respect to 9: 

R"{9) + G(R'{9),R(9),g ab , g ab , p , 5aM ) = 0. (62) 

G is a rather lengthy function of its arguments, non-linear in R and R'; for brevity we do not display it 
explicitly. All of the metric functions and their gradients appearing in (|62|l are evaluated along a given curve 
of integration, and hence are implicitly functions of 9. 9 ranges from to tt, and regularity of the surface 
F = about the axis requires R'(Q) — R'(tt) — 0. Integration of 1|62[1 therefore proceeds by specifying R at 
9 = (for instance), and then "evolving" R until either 9 = tt, or R diverges at some value of 9 < tt. If an 
AH exists, and assuming zq is inside the AH, then the AH can be found by searching for the (locally) unique 4 
initial value R(0) = Ro such that integration of l|62() ends at R(tt) = R-^, with R'(tt) = and R^ finite. For 
R(0) slightly larger than Rq (outside the AH), the integration will end at 9 — tt with R'(tt) > 0, indicating 
an irregular point on the surface; similarly, for R(0) slightly smaller than Ro (inside the AH) the integration 
will end with R'(tt) < 0. Therefore, if we can find a reasonable bracket about the unknown Rq, we can use 
a bisection search to find Rq. Currently, we find a bracket to search by testing a set of initial points, equally 
spaced in z at intervals of 3Az. This seems to work well in most situations, and the search is reasonably fast. 

We use a second-order Runge-Kutta method to integrate equation (|62f> . The metric functions appearing in 
G are evaluated using bilinear interpolation along the curve. 



V. IMPLEMENTATION 



In this section we describe the numerical code that we have written to solve the equations listed in Section 
IIIII Some details are deferred to Appendix lAl 

We use a uniform grid of size N p points in p by N z points in z, with equal spacing Ap = Az = h in the p 
and z directions. The value of a function / at time level n and location (i, j) within the grid, corresponding 
to coordinate (p,z,t) = ((i— l)Ap, (j — 1)A Z + z m - m , nAt), is denoted by /"•. For the temporal discretization 
scale we use At = Xh, where A is the Courant factor, which for the type of differencing we employ, should be 
less than one for stability; typically we use A = 0.3. The hyperbolic equations (|34I38|) are discretized using 
a second-order accurate Crank-Nicholson type scheme, whereby we define two time levels, t and t + At, and 
obtain our finite difference stencils by expanding in Taylor series about t = t + At/2. This gives the following 
second-order accurate approximation to the first derivative of / with respect to time 



fij-fij 9f(p,z,t) 



At dt 



+ 0(At 2 ) (63) 

t=t+At/2 



Second-order accurate approximations to functions and spatial derivative operators at t = t + At/ 2 are 
obtained by averaging the corresponding quantity, Q, in time: 

G n ■ + O n+1 

Wl ' J ^ J = Q(t + At/2, p, z) + O(At') (64) 

Thus, after discretization of the evolution equations using (|63|l and (|64|l . function values are only referenced 
at times t and t + At, even though the stencils are centered at time t + At/2. Specific forms for all the finite 
difference stencils that we use can be found in Appendix lAl 

We add Kreiss-Oliger dissipation [2^ to the evolution of equations of $, n, a and O (in addition to ip during 
partially constrained evolution), as described in Appendix lAl To demonstrate that this is essential for the 
stability of our numerical scheme, we compare in Fig. nt ne evolution of a from simulations without and with 
dissipation, but otherwise identical. 

The elliptic equations l|3UI33|l are solved using Brandt's FAS multigrid (MG) algorithm |U|33|, described 
in some detail in Appendix There are no explicit time derivatives of functions in these equations, and 
we discretize them at a single time level n (i.e. we do not apply the Crank- Nicholson averaging scheme for 
the elliptics). We use either a fully constrained evolution, solving for a,(3 p ,/3 z and ip using the constraint 



4 In a general collapse scenario, multiple inner horizons could be present, which would also satisfy 1621 augmented with the 
conditions i?'(0) = R' (ir) = 0. We want the outermost of these surfaces. 
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FIG. 1: The metric variable a from two sample Brill wave evolutions, demonstrating the effectiveness of dissipation. 
Both simulations were run with identical parameters, except the simulation displayed on the left was run without 
dissipation, while that shown on the right was run with a dissipation parameter e = 0.5 (see Appendix [XJ. The 
simulation without dissipation crashed at around t = 10, while the run with dissipation was stopped after several light 
crossing times, without showing any signs of instability. (Note: the computational domain used for these evolutions 
was p = [0, 10], and z = [—10, 10]; in the plots above we only show a subset of this domain. Also, the grid-lines drawn 
are for visual aid only; the size of the mesh visible in each frame is 64 x 128.) 



equations and slicing condition, or a partially constrained evolution where instead of using the Hamiltonian 
constraint to update ip, we use the evolution equation l|36() . Partially constrained evolution has proven to 
be useful due to the occasional failure of the MG solver in the strong-field regime (i.e. close to black hole 
formation). Use of the evolution equation for ip (rather than the Hamiltonian constraint) circumvents this 
problem in many instances; however, in certain Brill-wave dominated spacetimes, free evolution of ip is not 
sufficient to ensure convergence of the MG process. We are currently working to make the MG solver more 
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robust in these situations. 

The code is written in a combination of RNPL (Rapid Numerical Prototyping Language (3(J) and Fortran 
77. The hyperbolic equations are implemented in RNPL, which employs a point-wise Newton-Gauss-Seidel 
iterative relaxation scheme to solve these equations, while the MG solver is implemented in Fortran (see 
Appendix 1X1 for more MG details). A pseudo-code description of the time-stepping algorithm used is as 
follows: 

a) As an initial guess to the solution at time t+dt, copy variables 
from t to t+dt 

b) repeat until (residual norm < tolerance) : 

1: perform 1 Newton-Gauss-Seidel relaxation sweep of the 

evolution equations, solving for the unknowns at time t+dt 
2: perform 1 MG vcycle on the set of elliptic equations at 
time t+dt 
end repeat 

For the residual norm used to terminate the iteration we use the the infinity norm of the residuals of all 
updated unknowns. 



VI. TESTS 



In this section, we describe some of the tests we have performed to check that we are solving the correct set 
of equations. The first test consists of checking the equations against those derived with a computer algebra 
system |3l). By inputting the metric and coordinate conditions, the computer derived equations can then be 
subtracted from our equations and simplified. By finding that the differences simplify to 0, we can conclude 
that two sets of equations agree. 

For diagnostic purposes and as tests of the equations and of their discretization, we compute several 
quantities during the numerical evolution. The first is the ADM mass 24] 



M 



ADM 



1 

16tt 



lim 



b;a 



H a a . b ) N b dA, 



(65) 



where the integral is evaluated on a flat 3-space, i.e. with metric ds 2 = dp 2 + dz 2 + p 2 d<fi 2 . The spatial 
3-metric H a b is that from our curved space solution, but has its indices raised and lowered with the flat 
metric. Integrating around the boundaries of our numerical grid, the normal vectors N are ±d/dz and d I dp. 
After some algebra, the ADM mass becomes 



M A 



DM 



pV 4 



dp 

dp 
1 

Tp 



l 



dz. 



(66) 



A second set of quantities we calculate are the i^-norms of the residuals of the evolution equations for the 
extrinsic curvature l|22|) . which we denote E(Kab)- Because we do not directly evolve individual components 
of the extrinsic curvature, these residuals will not be zero; however, they should converge to zero in the limits 
as the discretization scale h — ► 0, and the outer boundary positions p max , z m ax, ~z m - m — > oo. Note that we 
include these last conditions because it is only in the limit r — > oo that our outer boundary conditions are 
fully consistent with asymptotic flatness. 

The convergence properties of our code are measured by computing the convergence factor, Q Ul associated 
with a given variable, it, obtained on grids with resolution h, 2h and Ah via 



K/t - U 2 h\\2 
\\u2h - Uh\\i ' 



(67) 



In particular, for the case of 0(h 2 ) (second order) convergence, we expect Q u — > 4 as h — > 0. 
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FIG. 2: Tests of spherically symmetric scalar collapse. The sequence of frames are from the evolution of an initial 
pulse in $ of the form 1521 with A§ — 0.02, e$ = 1, 7?$ = 7.0, A$ = 1.0, (p$,z<j) = (0,0). Shown are the functions 
r$ and — r 2 dip /dr. The output of the explicitly spherically symmetric code (with 2 10 radial grid points) is shown with 
solid lines, while the output of the axisymmetric code (N p = N z /2 — 2 8 ) is shown with dashed (a (p = 0, z > 0) slice) 
and dotted lines (az = slice). 

The first set of tests we present here are comparisons of $ and ip from the evolution of spherically symmetric 
initial data to the corresponding functions computed by a ID spherically symmetric code, the details of which 
are presented in Appendix In general, the results from the two codes are in good agreement. A sample 
comparison is illustrated in Fig. [21 which shows the scalar field obtained with the ID code as well as two 
radial slices of the corresponding solution calculated using the axisymmetric code. Note, however, that we 
do not expect exact agreement in the limit h — ► for a fixed outer boundary location, as the "rectangular" 
boundaries of the axisymmetric code are, in general, incompatible with precise spherical symmetry. 

In the second series of tests, we examine evolutions of Brill waves and non-spherical scalar pulses. Figs. [2HHI 
show results from two typical initial data sets, each computed using two distinct outer boundary positions. 
Each figure plots (a) the ADM mass Madm, (b) the ^2-norm of the residual of the pp component of the 
evolution equation for the extrinsic curvature E(K pp ), and, (c) the convergence factor ofip, as functions of 
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FIG. 3: Tests of Brill collapse using a initial pulse profile for a(0, p, z) of the form 1521 1 with = —3.0, — 0, A s = 1, 
£ct = 1, and (p^, Zct) = (0, 0). The evolution shown here corresponds to four crossing times and p max = z ma x = 10. The 
top frame shows the calculated ADM mass Madm- As the resolution increases, so does the level of mass conservation 
at early times, before energy has reached the outer boundary. The middle frame shows E(K PP ), the i?2-norm of 
the residual of the pp component of the evolution equation for the extrinsic curvature 12211 . At early times, before 
energy reaches the outer boundary, the residual decreases as the resolution increases. The bottom frame shows the 
convergence factor computed for the field tp. 



time (the convergence factor Q for other functions exhibit similar behavior as , and so for brevity we do not 
show them) . Here, one expects to see an improvement of the results — namely trends toward mass conservation 
early on, a zero residual, and a convergence factor of 4 — in the limits h — > and (p max , z max ) — > oo. After 
energy has reached the outer boundary, and to a lesser extent before (as is evident in the scalar field example 
in Figs.[S]and we fail to get consistency with the evolution equation (|22|l as h — > 0, for fixed (/O m ax, -z max ). 
This is a measure of the inaccuracy of our outer boundary conditions I|47I50[1 ; though the trends suggest that 
we do achieve consistency in the limit (p m ax, z max ) — > oo. 

Finally, we show some results of a simulation of black hole formation from the collapse of a spherically 
symmetric distribution of scalar field energy. Again, by looking at spherically symmetric collapse we can 
compare with the ID code (obtaining the same level of agreement as seen in the example in Fig. 0). However, 
here we want to show the behavior of our coordinate system (in particular the maximal slicing) in the strong- 
field regime, which demonstrates the need to incorporate black hole excision techniques and/or adaptive mesh 
refinement (AMR) before attempting any serious investigation of physics with this code. Fig. [7| shows plots 
of the ADM mass estimate (|rj5|) . an estimate of the black hole mass M aroa = ^ A/16n, where A is the area of 
the apparent horizon, and the minimum value of the lapse as a function of time from the simulation. Fig. [S] 
shows the conformal factor ip at several times, in the central region of the grid. 

Maximal slicing is considered singularity avoiding, because as the singularity is approached in a collapse 
scenario, the lapse a tends to 0, as demonstrated in Fig. [3 This effectively freezes the evolution inside the 
black hole, though it causes a severe distortion in the t — const, slices as one moves away from the black hole. 
This particular coordinate pathology is evident in Fig. [S] Recall from the 2-metric (|25ll that tp 2 determines 
proper length scales in the p and z directions; thus the rapid growth with time of ip shown in Fig. IHlmeans that 
a given coordinate area represents increasing proper area. Furthermore, the increase in magnitude of tp in the 
strong-field regime (which happens even when black holes do not form, and in non-spherical scalar field and 
Brill wave evolution, though not to the same extent as shown in Fig. |HJ) implies that our effective numerical 
resolution decreases in those regions, as some feature of the solution with a given characteristic size will span 
less of the coordinate grid. Thus, in the end, even though maximal slicing may prevent us from reaching a 
physical singularity, the "grid-stretching" effect is just as disastrous for the numerical code, preventing any 
long-term simulation of black hole spacetimes. For these reasons we will add black hole excision techniques 
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FIG. 4: Tests of Brill collapse with the same initial data and grid resolutions as those in Fig. |3] above, but with 
Pmax = Zmax = 20, and for only two crossing times. Notice the improvement in the behavior of the residual and mass 
aspect when energy reaches the outer boundary, as compared to the p max = Z max = 10 case above. 



and AMR before exploring physics with this code; our efforts in this regard are well underway, and will be 
described elsewhere. 



VII. CONCLUSION 



We have described a (2 + 1) + 1 gravitational evolution model which evolves axisymmctric configurations of 
gravitational radiation and/or a scalar field. A thorough battery of tests confirms that the correct equations 
are being solved. In particular, we have provided evidence that the code is second-order convergent, consistent 
and conserves mass in the limit where the outer boundary position goes to infinity. 

The unigrid code described here is the first step towards our long-term goal of studying a range of interesting 
theoretical and astrophysical phenomena in axisymmetry. These include gravitational collapse of various 
matter sources and gravitational waves, the corresponding critical phenomena at the threshold of black hole 
formation, head-on black hole collisions and accretion disks. To this end, we need to include support for 
angular momentum and additional matter fields in the code, as well as to add additional computational and 
mathematical infrastructure — adaptive mesh refinement, black hole excision and the capability of running in 
parallel on a network of machines. All of these projects are under development, and results will be published 
as they become available. 

Another goal of this project is to provide a platform from which to develop computational technology for 
3D work. In particular, we see development of AMR in axisymmetry as a precursor to its incorporation 
in 3D calculations. Likewise, accurate and stable treatment of boundary conditions presents a continual 
challenge in numerical relativity, and it is possible that we can develop an effective treatment of boundaries 
in axisymmetry that will generalize to the 3D case. 
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FIG. 5: Tests of an oblate, scalar pulse evolution with A$ = 0.15, = 0, A$ = 3, e$ = 3, and (p$, z$) = (0, 0). The 
tests are shown for four crossing times with p max = z max = 10. The results are similar to the tests for the Brill wave 
shown in Fig. |3] above, though notice that the scale of the residual E(K PP ) is about 2 orders of magnitude smaller 
than that of the Brill wave case. 



NSF PHY-0139980 as well as the financial support of Southampton College. EWH also acknowledges the 
support of NSF grant PHY-0139782. The majority of the simulations described here were performed on 
UBC's vn cluster (supported by CFI and BCKDF), and the MACI cluster at the University of Calgary 
(supported by CFI and ASRA). 



APPENDIX A: THE ELLIPTIC SOLVER AND FINITE DIFFERENCE APPROXIMATIONS 

In this appendix we briefly mention some aspects of our multigrid (MG) routine, and list the set of finite 
difference approximations that we use. 

The constraint equations (|3C)133|I arc four elliptic equations which, for a fully constrained system, must be 
solved on every time slice (i.e. spatial hypersurface) . As such, it can be expected that the time taken by a 
given evolution will be dominated by the elliptic solver and hence we look for the fastest possible solver. 

Currently, multigrid methods are among the most efficient elliptic solvers available, and here we have 
implemented a standard Full Approximation Storage (FAS) multigrid method with V-cycling (see |U|3^|) 
to solve the four nonlinear equations simultaneously. (When using the evolution equation for ip in a partially 
constrained evolution, we use the same multigrid routine described here, except we only solve for the three 
quantities a,(3 p and (3 Z during the V-cycle; if) is then simply considered another "source function".) 

A key component of the MG solver is the relaxation routine that is designed to smooth the residuals 
associated with the discretized elliptic equations. We use point-wise Newton-Gauss-Seidel relaxation with 
red-black ordering (see Fig. 0), simultaneously updating all four quantities (a, /? p , {3 Z , ip) at each grid-point 
during a relaxation sweep. In addition to its use for the standard pre-coarse-grid-correction (pre-CGC) and 
post-CGC smoothing sweeps, the relaxation routine is also used to solve problems on the coarsest grid. We 
use half-weighted restriction to transfer fields from fine to coarse grids and (bi)linear interpolation for coarse 
to fine transfers. We generally use 3 pre-CGC and 3 post-CGC sweeps per V-cyc\e, and likewise normally 
use a single U-cycle per Crank-Nicholson iteration. 

One complicating factor here is the treatment of the boundary conditions Eqs. (|39I42|) and Eqs. H47I50JI ■ In 
accordance with general multi-grid practice, we view the boundary conditions as logically and operationally 
distinct from the interior equation equation. The outer boundary conditions can generally be expressed as 



r X rs constant 



(Al) 
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FIG. 6: Tests of an oblate, scalar pulse evolution with the same initial data and grid resolutions as those in Fig. |S] 
above, but with p max = w = 20, and for two crossing times. Again, as with the Brill wave example, it is evident 
that there are two factors that contribute to a non-zero residual E(K PP ) — the closeness of the outer boundary, and 
the discretization scale h. 




FIG. 7: The ADM mass estimate Madm of the spacetime and area-mass estimate M area = y/ A/16 of the black 
hole (left), where A is the area of the apparent horizon, and the natural logarithm of the minimum value of lapse 
(right) during a scalar field collapse simulation. The initial scalar field profile, $(0,p, z) is of the form 1521 with 
A$ = 0.35, R<s> = 0, A$ = 1, eg, = 1 and (p$, z$) = (0, 0). The outer boundary is at p max = z max = -z m in = 10, and 
the size of the numerical grid is 256 x 512. An apparent horizon was first detected at t w 4, hence the M aroa curve only 
starts then. At intermediate times (between roughly t = 2 and t — 7) we see an exponential "collapse" of the lapse 
(the minimum of which is at (p, 2) = (0, 0)); at later times this behavior ceases in the simulation, as a consequence of 
increasingly poor resolution in the vicinity of the black hole caused by "grid-stretching" — see Fig. |H] below. This also 
adversely affects the accuracy of the area-mass estimate (it actually begins to decrease at late times) as the coordinate 
region occupied by the AH shrinks. 



where X e {1 — a, 1 — if>, f3 z , f3 p }. Taking the derivative of IjAlfl with respect to r, we arrive at the differential 
form that is applied on the outer boundaries of the computational domain {p = p max , z — z max , z = —Zmax): 

X - pA p — zl z = 0. (A2) 

On the z-axis {p = 0), the conditions Eqs. (|39I42|I take one of the following two forms: 



A i!P = 



(A3) 
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FIG. 8: The conformal factor ijj near the origin, at four different times from the simulation described in Fig. [7] The 
height of the surfaces represent the magnitude of ip, and the scale of each image is the same. The smallest value of 
ip shown in each frame is « 1.5, while at t — 15, ij) reaches a maximum of « 5.8 at the origin. The 2-metric has the 
form ip 4 (dp 2 + dz 2 ), hence the larger ip is, the larger the physical area represented by a given coordinate cell (the lines 
drawn here are coincident with the actual grid-lines of the simulation). This implies that the effective resolution of 
a given coordinate patch decreases as if) increases. The net result is that as gravitational collapse proceeds, central 
features of the solution become very poorly resolved within the grid, adversely affecting the accuracy of the solution. 
This is quite evident at t = 15, where noticeable asymmetries have developed in tp (recall that this is collapse from 
spherically symmetric initial data). 



or 

At = (A4) 

Equations of the former form are discretized using an 0{h 2 ) backwards difference approximation to the p 
derivative. 

The interior and boundary differential equations are solved in tandem via the multigrid approach: 

1. The residual is smoothed using some number of relaxation sweeps. For the interior, Eqs. Q3UI32I) are 
relaxed using red-black ordering as discussed above. After each call of this relaxation routine, a second 
routine that "relaxes" the boundary points is called |35j . 

2. For quantities restricted from a fine to a coarse grid, discrete forms of i|A3(l and l|A4(l are applied during 
the V-cycle. At the other boundaries, straightforward injection is used. 
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FIG. 9: Diagram illustrating the order in which relaxation occurs. Shown here is a 5 x 9 grid spanning p x z with the 
grid points represented by filled circles. The Xs denote "red" (interior) points which are visited first. The Os denote 
"black" (interior) points which are visited next. Finally, the triangles denote boundary points which are visited last. 

The key idea here is to ensure that the boundary relaxation process does not substantially impact the 
smoothness of the interior residuals, because it is only for smooth residuals that a coarsened version of a 
fine-grid problem can sensibly be posed. 

Finally, in Table U we show all of the difference operators we use to convert the differential equations listed 
in IIIII to finite difference form, using the Crank-Nicholson scheme described in Section El In addition, as 
discussed in Section E] we use Kreiss-Oliger dissipation to maintain smoothness in the evolved fields. 
Specifically, we add the Kreiss-Oliger filter to discretized evolution equations 

A t A n = n t f n (...) (A5) 
by replacing the Crank-Nicholson time difference operator A t with A^: 

A\A n = ix t p (. . .) (A6) 
Empirically, we find that a value of e = 0.5 generally keeps our fields acceptably smooth. 



APPENDIX B: THE SPHERICALLY SYMMETRIC MODEL 



One simple test of the code compares the results for spherically symmetric initial data with the output of 
a code which explicitly assumes spherical symmetry. Here we present the equations for this ID code. The 
spacetime metric is: 

ds 2 = -(a 2 + ij/f3 2 )dt 2 + 2ip A /3dtdr + A (dr 2 + r 2 dft 2 ) , (Bl) 

where a, (3 and ip> are functions of r and t, dft 2 is the line element on the unit 2-sphere, and (3 is the radial 
component of the shift vector (i.e. (3 l = (/3,0, 0)). Adopting maximal slicing to facilitate direct comparison 
to the axisymmetric code, we have 

K l j = diag(iT r (r, t), 0, 0). (B2) 
Then a sufficient set of equations for the coupled Einstein-massless-scalar system is [34| 

V/' + ^ + 2n [4> 2 + n 2 ] i> + ^ (K) 2 A = (B3) 
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A x Ai 

AlAi 
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A x (Ai/x) 



(A i+ i - Ai-x) I (2ft) 
(-3A, + 4A I+ i - A l+2 ) I (2ft) 
(3Ai - 4Ai_! + Ai- 2 ) I (2ft) 
(Ai+i -2^ + ^-1) /(ft 2 ) 

(Ai+i - A 8 _i)/(:r 2 +1 -x 2 _ x ) 
2 (xi-iAi+i - Xi+iAi-i) I (x 2 +1 
A^A^O/x] 16 [sCi_i/ 2 A+i - 2xiAi + ac i+ i/ a j4.i_i] /(x 2 +1 - x 2 ^) 
HtA n (A n+1 + A n ) /2 

A t A n (A n+1 - A n ) I (Aft) 

A\A n [A t + eh 3 / (16A) A xxxx ] A™ 
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L + o(ft 2 ) 
l + o(h 2 ) 
l + o(h 2 ) 
*l + o(h 2 ) 
xxx \ l + o(h 2 ) 
*\> + o(h 2 ) 



{A/x), x \ i + 0(h 2 ) 
(A, x /x), x \ t + 0(h 2 ) 
A\ n+1/2 +0(\ 2 h 2 ) 
A, t \ n+1/2 + 0(X 2 h 2 ) 
P(\ 2 h 2 ) 



A. 
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TABLE I: Finite difference operators and their correspondence to differential operators. Here, A™ is an arbitrary grid 
function defined vis, j4j = -A(Xmin ~\~ (t 1)^; ^min + (n — l)(Aft)), where ft and Aft are the spatial and temporal grid 
spacings, respectively, x denotes either of the two spatial coordinates p or z, with the dependence of A™ on the other 
suppressed. The parameter e represents a user-specifiable "amount" of Kreiss-Oliger dissipation. 
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(B6) 
(B7) 

(B8) 
(B9) 



Here dot and prime denote derivatives with respect to t and r, respectively. 

The evolution equations l|B7IB9(l are discretized using an 0(h 2 ) Crank-Nicholson scheme. Eqs. (£2; 
and l|B4|) are similarly discretized using 0(h 2 ) finite difference approximations, then solved iteratively for 
ifj and K r r at each time step. Once $, n, ip and K r r have been determined, a and (3 are found from 0(h 2 ) 
finite-difference versions of Eqs. 1B5J1 and l|B6|l . The code is stable and second-order convergent. 



[1] L. L. Smarr, Ph.D. dissertation, University of Texas at Austin, unpublished (1975) 
[2] K. R. Eppley, Ph.D. dissertation, Princeton University, unpublished (1977) 

[3] S. L. Shapiro and S. A. Teukolsky, "Collisions of Relativistic Clusters and the Formation of Black Holes", Phys. 
Rev. D45, 2739 (1992) 

[4] P. Anninos, D. Hobill, E. Seidel, L. Smarr and W. Suen, "The Collision of two black holes," Phys. Rev. Lett. 71, 
2851 (1993) gr-qc/9309016 . 

[5] J. Baker, A. Abrahams, P. Anninos, S. Brandt, R. Price, J. Pullin and E. Seidel, "The Collision of Boosted Black 
Holes", Phys. Rev. D55, 829 (1997) 



19 



[6] P. Anninos and S. Brandt "Head-on Collision of Two Unequal Mass Black Holes," Phys. Rev. Lett. 81, 508 (1998) 
gr-qc/9806031 . 

[7] T. Nakamura, "General Relativistic Collapse of Axially Symmetric Stars Leading to the Formation of Rotating 

Black Holes," Prog, of Theor. Physics 65, 1876-1890 (1981). 
[8] R. F. Stark and T. Piran, "Gravitational Wave Emission From Rotating Gravitational Collapse," Phys. Rev. Lett. 

55, 891 (1985). 

[9] T. Nakamura, K. Oohara and Y. Kojima, "General Relativistic Collapse of Axially Symmetric Stars", Prog. 

Theor. Phys. Suppl. 90, 13 (1987) 
[10] M. Shibata, "Axisymmetric Simulations of Rotating Stellar Collapse in Full General Relativity — Criteria for 

Prompt Collapse to Black Holes" Prog. Theor. Phys. 104, 325 (2000). 
[11] A. M. Abrahams, G. B. Cook, S. L. Shapiro and S. A. Teukolsky, "Solving Einstein's equations for rotating 

space-times: Evolution of relativistic star clusters," Phys. Rev. D 49, 5153 (1994). 
[12] S.L. Shapiro and S.A. Teukolsky, "Formation of Naked Singularities: The Violation of Cosmic Censorship," Phys. 

Rev. Lett. 66, 994 (1991). 

[13] D. Garfinkle and G. C. Duncan, "Numerical evolution of Brill waves." Phys. Rev. D 63, 044011 (2001) 
gr-qc/0006073 . 

[14] M. Alcubierre,S. Brandt, B. Brugmann,D. Holz,E. Seidel,R. Takahashi,J. Thornburg "Symmetry Without Sym- 
metry: Numerical Simulation Of Axisymmetric Systems Using Cartesian Grids " Int. J. Mod. Phys. D 10, 273 
(2001). 

[15] S. Brandt, J. A. Font, J.M. Ibanez, J. Masso and E. Seidel, "Numerical Evolution of Matter in Dynamical Ax- 
isymmetric Black Hole Spacetimes", Comput. Phys. Commun. 124 169 (2000) 

[16] S. Brandt and E. Seidel, "The Evolution of Distorted Rotating Black Holes I: Methods and Tests", Phys. Rev. 
D52, 856 (1995) 

[17] S. Brandt and E. Seidel, "The Evolution of Distorted Rotating Black Holes H: Dynamics and Analysis", Phys. 
Rev. D52, 870 (1995) 

[18] S. Brandt and E. Seidel, "The Evolution of Distorted Rotating Black Holes III: Initial Data", Phys. Rev. D54, 
1403 (1996) 

[19] A. M. Abrahams and C. R. Evans, "Trapping A Geon: Black Hole Formation By An Imploding Gravitational 
Wave.," Phys. Rev. D. 46, 4117 (1992). 

[20] A. M. Abrahams and C. R. Evans, "Critical behavior and scaling in vacuum axisymmetric gravitational collapse," 
Phys. Rev. Lett. 70, 2980 (1993). 

[21] K. Maeda,M. Sasaki, T. Nakamura,S. Miyama "A New Formalism of the Einstein Equations for Relativistic Ro- 
tating Systems." Prog. Theor. Phys. 63, 719 (1980). 

[22] K. Macda, "[(2+l) + l]-Dimcnsional Representation of the Einstein Equations," in Proceedings of Third Marcel 
Grossmann Meeting on General Relativity, ed. Hu Ning, (Science Press, 1983). 

[23] R. Geroch, "A Method For Generating Solutions Of Einstein's Equations," J. Math. Phys. 12, 918 (1971). 

[24] C.W. Misner, K.S. Thorne and J.A. Wheeler, Gravitation, New York, W.H. Freeman and Company (1973) 

[25] J.M. Bardeen and T. Piran, "General Relativistic Axisymmetric Rotating Systems: Coordinates and Equations", 
Phys. Rep. 96 205 (1983). 

[26] T. Appelquist, A. Chodos, and P.G.O. Freund, "Modern Kaluza-Klein Theories," Addison- Wesley, Menlo Park, 
1987). 

[27] J. Thornburg, "Finding Apparent Horizons in Numerical Relativity", Phys. Rev. D54, 4899 (1996) 

[28] M. Alcubierre, S. Brandt, B. Brugmann, C. Gundlach, J. Masso, E. Seidel and P. Walker, "Test Beds and 

Applications for Apparent Horizon Finders in Numerical Relativity", Class. Quant. Grav. 17 2159 (2000) 
[29] Kreiss, H.-O., and Oliger, J., "Methods for the Approximate Solution of Time Dependent Problems", Global 
Atmospheric Research Program Publication No. 10, World Meteorological Organization, Case Postale No. 1, 

CH-1211 Geneva 20, S witzerland (1973). 

[30] Software available from http://laplace.physics.ubc.ca/Members/matt/Rnpl/index.html 
[31] We use Waterloo Maple along with a tensor package written by one of the authors (MWC). 

[32] A. Brandt, "MultiLevel Adaptive Solutions to Boundary Value Problems," Math, of Computation 31, 333-390 
(1977). 

[33] A. Brandt, "Guide to Multigrid Development," in Lecture Notes in Mathematics 960, 220-312, (Springer- Verlag, 
New York, 1982). 

[34] MW. Choptuik, "A Study of Numerical Techniques for Radiative Problems in General Relativity," Ph.D. thesis, 
The University of British Columbia, unpublished (1982). 

[35] A subtle point here concerns the "relaxation" occurring on the boundaries. The finite-difference equations on the 
boundary yield algebraic equations which determine the given fields there "exactly." The subtlety arises because 
these algebraic conditions couple neighboring boundary points and thus despite the fact that we solve these 
equations exactly, the residual (as computed after the entire boundary is "relaxed") will not be identically zero. 
We find this procedure suffices to keep residuals sufficiently smooth over both interior and boundary domains. 



